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Abstract 

We develop an efficient method to calculate the third-order corrections to 
the self-energy of the hole-doped two-dimensional Hubbard model in space- 
time representation. Using the Dyson equation we evaluate the renormalized 
spectral function in various parts of the Brillouin zone and find significant 
modifications with respect to the second-order theory even for rather small 
values of the coupling constant U. The spectral function becomes unphysical 
for U ~ W, where W is the half-width of the conduction band. Close to the 
Fermi surface and for U < W, the single-particle spectral weight is reduced 
in a finite energy interval around the Fermi energy. The increase of U opens 
a gap between the occupied and unoccupied parts of the spectral function. 

Introduction 

The hole-doped two-dimensional (2-D) Hubbard model continues to attract consid- 
erable attention as it is believed to capture the essential features of "underdoped 
cuprates" |l|]. In particular, the spectral properties of the model at small doping 
and the coupling strength U > 2W, where W is the half-width of the conduction 
band, are thought to be relevant for the angular resolved photoemission spectroscopy 
(arpes) experiments. The arpes data are used to study the frequency and mo- 
mentum dependence of the single-particle excitations of 2-D electrons in C11O2 lay- 
ers @, U but the interpretation of the experimental results is quite difficult and 
often controversial. On the one hand, the observed spectral features do not fit any 
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simple conceptual framework; on the other hand, reliable theoretical results are not 
at hand. Thus, it is of interest to study in some detail the effect of correlation on 
the spectral function of the hole-doped 2-D Hubbard model. 

The weak-coupling analysis of renormalized single-particle excitations has been 
presented in a number of papers P}|-||13||, which treat U as an expansion parameter 
and consider the effects of correlation, doping, and temperature in various parts 
of the Brillouin zone. These papers show that correlation changes significantly the 
single-particle spectral properties even for relatively small values of U, and the 
results exhibit a number of interesting features which are also seen in cuprates. 
However, a quantitative comparison with the experimental data is difficult since 
most weak-coupling theories become unreliable for U > W. 

The breakdown of the weak-coupling schemes based on truncated perturbation 
expansions is not immediately discernible from the spectral function, but is signalled 
for U W by the negative compressibility |TJ| and the rapid deviation of the Fermi 



volume vp from the values required by the Luttinger theorem [14]]. The theories 
based on an infinite summation of selected classes of diagrams |fT5| , [16], [T7J are also 
unreliable for large values of U, because they overemphasize the spectral weight of 
quasiparticle peaks and do not produce the Hubbard side-bands, which are typical 
of strong local correlations. Thus, the perturbational results obtained so far give 
some insights into the properties of the Hubbard model but do not allow a consistent 
description of correlated electrons at intermediate or large values of U which one 
needs to discuss the cuprates or make a comparison with the t-J model. 

The weak-coupling approach to the Hubbard model poses a number of ques- 
tions which should be considered if the results are to be extrapolated into the large-t/ 
regime. What is the range of validity of the asymptotic U expansion? What is the 
importance of the terms that are neglected in the self-energy expansion? Is it pos- 
sible to use a finite-order perturbation theory for the values of U such that the 
low-energy excitations of the Hubbard model and the t-J model look similar? Are 
the weak-coupling results obtained by perturbative methods representative of the 
strong-coupling limit? And finally, does the low-energy physics of the 2-D Hubbard 
model, as defined by some weak-coupling scheme, produce the right phenomenology 
for underdoped cuprates? 

To answer these questions one would have to examine the general structure of 
the perturbation expansion or compare the weak-coupling solution with the exact 
one, which is not possible at present. Some insight, however, can be obtained 
by extending the perturbation expansion beyond the 2nd order and studying the 
stability of a truncated series with respect to higher-order corrections. Here, we 
calculate the momentum-dependent single-particle self-energy up to the 3rd order, 
and show that the individual 3rd-order diagrams are large but that the total 3rd- 
order contribution is much smaller than the 2nd-order one. This approximate mutual 
cancellation of 3rd-order diagrams holds for any doping and not just for zero doping, 
where it is a consequence of electron-hole symmetry. 

Once the 3rd-order contribution becomes comparable to the 2nd-order one, 
which happens here for U ~ W, the truncated perturbation series leads to un- 
physical results. In the "physical" range, U < W, we find a number of interesting 
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features which offer additional insight into the anomalies of the 2-D hole-doped 
Hubbard model. However, the error of neglecting the higher-order terms might 
become significant for larger values of U even within the physical range [|11|. It 



would be interesting to see whether the 4th order contribution improves the pertur- 
bation theory and extends it to experimentally relevant values of U, or whether it 
renders the weak-coupling approach useless. The technical problems involved with 
such calculations are not much different from the problems encountered in the 3rd- 



order calculations |19| . The Matsubara summations for the 4th-order diagrams are 
straightforward, if tedious, and the 4th-order momentum summations follow from 
the same numerical strategies which solve the 3rd order. Thus, the 3rd-order calcu- 
lations can be considered as a small but necessary step in our efforts to clarify the 
properties of the Hubbard model up to intermediate values of the coupling strength. 

We should also remark that the accurate characterization of the weak-coupling 
regime might be of some interest for the approximate schemes which interpolate 
between the small-?/ and the large-?/ limit of the model. In the case of the infinite- 
dimensional Hubbard model and the single-impurity Anderson model, interpolations 
like that [20|— [23] come very close to the exact solution. However, on a 2-D lattice 
one deals with the anisotropic k space, and the interpolation schemes might be 
difficult to construct. 

The present paper is organized as follows. First the retarded 2nd- and 3rd- 
order self-energies are expressed in terms of multiple momentum integrals. Then we 
introduce the space-time representation and use the Fast Fourier Transform (fft) 
algorithm to evaluate the self-energy as a function of energy for all points in the 
Brillouin zone. The relative importance of the 2nd- and 3rd-order terms obtained 
by fft is analyzed and the stability of the weak-coupling solution is discussed. 
Next, the spectral properties of the model are calculated for low temperatures, for 
a given value of the chemical potential, and for various points in the Brillouin zone. 
Finally, the spectral features and their relevance for the experimental data are briefly 
discussed. The calculations are explained in detail in the Appendix. 



Calculations 

To start with, the Hubbard Hamiltonian is written as 

H = } J tijCj .Cj <T - yL}n i(T + U yjn^- «n iT , (1) 

i,j,cr i,tr i=l 

where Uj is the nearest-neighbor hopping integral, q ct (ctj destroys (creates) an elec- 
tron at site Rj with spin a, rii a = c\ a Ci a is the electron number operator, U is the 
local electron-electron interaction, and denotes the ensemble average over 

the eigenstates of the full Hamiltonian (p. The parameter // = \i — «rt_„ y> U 
is the "effective chemical potential", \i being the chemical potential proper. The 
energy of an unrenormalized single-particle excitation propagating with wave vector 
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k, counted from the Fermi level, is a;£ = — fi f , where ek = —2[t x cos(k x a x ) + 
tyCos(k y a y )}. The Fermi momentum is denoted by kp. We consider the fluctua- 
tions above a mean-field-like paramagnetic state, in which the number of particles 
coincides with the exact particle number, <<< ^> — « >>> = n e /2, and assume 
t x = ty = t, the half-bandwidth being W = 4t. 

The self-energy diagrams are generated by the usual Matsubara imaginary- 
time perturbation expansion with respect to the last term in Eq. (jl]). All the 
2nd- and 3rd-order diagrams are shown in Fig. 1, where the dashed line represents 
the local interaction (-U) and the full line stands for the unperturbed propagator 
G^iiuJn) = {iuj n — u}^)~ 1 , defined with the first two terms in Eq. (jl]). Note that there 
are no self-energy diagrams that are k and u independent, the so-called one-legged 
diagrams. Because of the separation of the Hamiltonian as in ([!]), all one-legged 
diagrams of all orders add up to zero. Equivalently one could say that they have 
been hidden in the "effective chemical potential" y! = fi — <<C n^ a >>> U . 

The self-energy calculations for a finite lattice with periodic boundary condi- 
tions and a discretized time axis simplify considerably in the space-time represen- 
tation m . The 2nd- and the 3rd-order retarded proper self-energy contributions are 
given by the expressions 

Eg } (t) = -iU 2 Q{t) [a 2 (R,£)b*(R,t) + b 2 (R, t) a*(R, t)] , (2) 
Eg } (t) = U 3 Q(t) {[a(R,t)wi(R,t) -b(R,t)w*(R,t)] 

- [a*(R,t)w 2 (R,t)+b*(R,t)w 3 (R,t)]} , (3) 

where a(R, t) and b(R, t) are the standard |MJ double-time Green's functions of H , 



a(R, t) = < 4(0) pa(f) > = 1 e i(k R -^ f , (4) 

k 

b(R, t) = < CB.it) 4(0)> = ^E e i(k - R ~ <t] f(-^k) , (5) 

k 

while the functions Wj(R, t) are space-time convolutions of functions composed of 
products of a(R, i)'s and b(R, t)'s, given by expression (|24"D in the Appendix. Instead 
of evaluating these convolutions directly, we decouple them by a pair of Fourier 
transforms, as shown in (p5|). (The time variable t should not be confused with the 
hopping amplitude t.) 

The retarded self-energy in the energy-momentum representation, £k(o>) = 
£k(cc + ^0 + ), is then given by the inverse Fourier transform 

/■oo 

£ k (o;) = J>- ik - R / dte^ Sg J (0 + Sg?(0 • (6) 
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The self-energy calculation is thus reduced to a sequence of Fourier transforms, which 
can be evaluated very efficiently by the fft technique. In this paper we consider a 
lattice with N g = 256 x 256 sites and define discretized momenta in the quadratic 
Brillouin zone as k = (k x , k y ), where k x , y = Ak(l x>y — 1) with Ak = 2ir/ ^/N^ and 
lx,y = !)•••> V^fl" ^he T P°i n t is at k = (0, 0), the X point at k = (n, 0), the M 
point at k = (tt/2, 7r/2), and the Z point at k = (n, 7r). 

Having found Ek(u;,T), we calculate the spectral function A^(uj,T) from the 
Dyson equation, 

A k (a;,T) = -ilm— — : \ — , , , (7) 

7r + i0 + — a;^ — L k (c<j, T) 

and obtain the renormalized density of states, 

p(^ T )-^J2M^n (8) 

k 

and the renormalized particle number, 

/oo 
duj{(uj) p(u, T) . (9) 
-oo 

Eq. (0), together with Eqs. @ and (0), establishes the functional dependence 
of the renormalized particle number n e on the "effective chemical potential" fi' , with 
U and T as parameters, everything scaled in units of W: 

n e = N^'/W- U/W, k B T/W) (10) 
On the other hand, the pure chemical potential [i is given by 

, u 

[L = p. + —n e 

= fi' + jAr(ji'/W;U/W,k B T/W) (11) 



Eqs. (pTOf ) and flllD establish the dependence of n e on fx or vice versa, with both 
quantities given parametrically as functions of y! . In this way one does not have 
to really solve Eq. @ as a self-consistency equation for n e (/i), it suffices to just 
evaluate Eqs. ( |10) and (|TT| ) for a number of close-lying values of \i! . 



One should note that the n-consistency is here forced upon the approximate 
proper self-energy, H^\uj) + Ilj^u;), and can therefore be attained only approxi- 
mately. 
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The peaks of A k (u) give the dispersion of quasiparticle excitations, u^; the 
renormalized Fermi surface is defined by the set of points in the momentum space 
at which = 0. (At the Fermi surface, A^(uj) has a singularity at the Fermi 
energy Ep.) According to the Luttinger theorem, the number of k points enclosed 
by the Fermi surface (i.e. the Fermi volume vf) should coincide with n e . In our 
approximate treatment, which is based on the 3rd-order self-energy and the Dyson 
equation, we find vf > n e , but the relative difference between vf(U,{i) and n e (U, /i) 
for U < W is very small. 

In what follows, we first discuss £ k (co>,T) and then the ensuing spectral func- 
tion, the density of states, the renormalized dispersion, and the Fermi surface of 
the model for the temperature k B T = t/250, U = 3.5t, and // = 0.2t, which cor- 
responds to n e = 0.96. The same numerical program, available from the authors 
upon request, returns the self-energy in the irreducible wedge of the Brillouin zone 
for any other value of T, U, and \x. 

Results and discussion 

In Figs. 2 and 3 we show the real and the imaginary part of £^(u;), respectively, 
plotted versus to for several momenta along the r-X-Z-M-T cuts through the Bril- 
louin zone. The dotted and the dashed-dotted lines show the two 3rd-order contri- 
butions, Ep^(u;) and £p P (u;), respectively, while the full line is their sum, T,^\u>). 
The magnitudes of the individual 3rd-order contributions are about the same as the 
magnitude of H^\u) (see below and in Ref. ||) but the total 3rd-order contribution 
is relatively small, except in a narrow frequency range where the approximate can- 
cellation of electron-electron and electron-hole terms does not occur. This behavior 
indicates that the correct solution of the hole-doped model requires all 3rd-order 
diagrams, and can not be obtained by the partial summation of electron-electron or 
electron-hole diagrams. 

The functional form of £^ (a;) is very much k-dependent. Far off the Fermi 
surface ImS^ (u)) does not show much structure in the low-energy region and its 
slope around Ep is always small (see Figs. 3 and 5 for data corresponding to T 
and Z points). On the other hand, for k ~ kp, ImE^ (u;) acquires a pronounced 
minimum at Ep (see Figs. 3 and 5 for data corresponding to X and M points). In 
the low-energy region, the ImE^ 3 ' (w) and ImE^(w) have opposite sign, and both 
vanish at Ep [^5| with a zero slope. 

The properties of the "full" self-energy, t!^\uo) + Tj^\uj), are summarized 
in Fig. 4, where ImSk(cu) is plotted versus uj for the same k's as in Fig. 3. The 
low-energy behavior at T, X, Z, and M is displayed in Fig. 5. For comparison, 
the individual 2nd- and 3rd-order contributions are shown as well. The low-energy 
features of the Re£k(u;) are n °t changed much by the 3rd-order renormalization, re- 
gardless of k. On the other hand, the low-energy properties of ImSk(u;) are strongly 
influenced by the 3rd-order corrections. For example, close to X and M points, 
the 2nd- and the 3rd-order terms nearly cancel and make the low-energy part of 
ImEk(co') rather flat in an extended interval around Ep. The relative contribution of 
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the 2nd- and the 3rd-order terms depends on the coupling strength and for a given 
\x there is always a critical correlation U c (fi, k c ) such that the coefficient of the u 2 
term in ImEk c (o;) vanishes at some k c . For U < U c the curvature of ImEi c (o;) at Ep 
is negative for all k and the ensuing spectral weight is always positive. However, 
for U > U c the curvature of ImSk(cu) at Ep becomes positive at some points in 
the Brillouin zone, and the corresponding spectral weight becomes negative. This 
behavior indicates that the 3rd-order U expansion breaks down for U ~ W, and 
that E k (u) provides an accurate renormalization for small values of U only. 

For U close to but less than U c , the curvature of ImS k (w) at E F is found to 
be the smallest for k ~ kp. Thus, in the presence of correlations the Fermi surface 
properties assume qualitatively new features due to the u 3 and higher-order self- 
energy terms. Unfortunately, these non-Fermi-liquid (nfl) features require large 
values of U and can not be properly discussed without the higher-order terms in the 



expansion. Leaving the 4th-order renormalization for future studies [19j , we consider 
here the 3rd-order renormalization and discuss spectral properties for U < U c ~ W. 

The variation of ^4k(^) along r — > X, X —>■ Z and Z — > T cuts in the Brillouin 
zone is shown in Fig. 6. The 3rd-order spectral function, like the 2nd-order one |J, 
assumes at all wave vectors a typical shape with a low-energy quasiparticle peak 
and a high-energy incoherent background. The correlation reduces the quasiparticle 
spectral weight and enhances the incoherent background of Ak(u). However, for 
U = 3.5t and n e = 0.96 the transfer of spectral weight out of the low-energy region 
is small, and the quasiparticle peak is not fully separated from the incoherent back- 
ground, except for k's which are far off the Fermi surface. At the Fermi surface, the 
singular quasiparticle peak can be represented as A kp (uj) = Z^ F 5{uj — E F ), where 
Zk describes the reduction of the quasiparticle weight due to self-energy corrections. 
For k ^ hp the quasiparticle peak broadens, becomes asymmetric, and attains a 
maximum at some finite frequency. The asymmetry of A-^{u) is caused by an in- 
coherent background, which slows down the high-a; decay, and the NFL behavior at 
small uj. A new feature due to the 3rd-order, which begins to emerge for U ~ W, is 
the suppression of spectral weight in a finite energy interval around Ep. For k ~ hp, 
the singular quasiparticle peak at the Fermi energy is separated from the occupied 
and unoccupied incoherent states by a pseudogap. For k < kp a real gap appears 
between the (occupied) quasiparticle peak and the (unoccupied) incoherent states. 
(The small background seen in the numerical data is due to the damping factor used 
for the real- axis propagators.) 

The renormalized quasiparticle dispersion, u^, is obtained by tracing the max- 
imum of v4k(co>) across several cuts in the momentum space and is shown in Fig. 
7. Circles show and full line gives the unrenormalized dispersion, u;£, for the 
same number of particles. Correlation reduces the overall width of the quasiparticle 
band and extends the flat dispersion around (n, 0). Essentially the same dispersion 
is obtained by solving the secular equation 

^k(T) - (e k - n') - XkMT), T] = . (12) 
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The inspection of Figs. 6 or 7 shows that the M-point quasiparticle peak is somewhat 
above Ep, while the X-point peak is right at Ep. Thus, the renormalized Fermi 
surface, defined by the set of k points such that the quasiparticle peak is at E F , has 
a different topology than the non-interacting Fermi surface with the same number 
of holes. The renormalized Fermi surface resembles the tight-binding result which 
adds to tij the next-nearest-neighbor hopping t^. 

The renormalized density of states p(u) calculated for U = 3.5t and n e = 0.96 
is plotted in Fig. 8 as a function of u. The transfer of the low-energy spectral weight 
out of the low-cu region is clearly seen but U < W is not sufficient to generate the 
incoherent Hubbard side-bands. The density of states at E F is not much enhanced 
from the unperturbed value despite the reduced dispersion, because the quasiparticle 
weight is reduced by Z^. One has a metal but a strange one. 

In summary, we developed an efficient method to perform the momentum sum- 
mations in the self-energy expansion for the hole-doped 2-D Hubbard model and 
evaluated the coefficients of the U 2 and the U 3 terms. The third-order corrections 
lead to additional anisotropy and asymmetry of the spectral function and give rise to 
new features with respect to the 2nd-order renormalization. For large enough U we 
find strong reduction of the single-particle spectral weight around Ep. The quasi- 
particle dispersion is also reduced and the saddle point singularities are extended. 
The Fermi surface obtained by the 3rd-order renormalization for a hole-doped sys- 
tem is closed around the Z point. At the Fermi surface, the quasiparticle peak and 
the incoherent background are separated by a pseudogap. Off the Fermi surface we 
find a finite region of negligible spectral weight between the quasiparticle peak and 
the Fermi energy. This behavior hints to a possible scenario in which the system is 
metallic with a narrow Kondo-like resonance at the Fermi energy, and a pseudogap 
for incoherent excitations. However, the anomalous spectral features which we ob- 
tain here are due to the cancellation of the 2nd- and the 3rd-order self-energy terms, 
and the stability of this result with respect to the 4th-order renormalization remains 
to be seen. The 3rd-order perturbation expansion shows that the solution of the 
hole-doped model can not be obtained by partial summation of diagrams, and that 
the properties inferred from the 2nd-order theory seem to be qualitatively correct 
for U < W only. 



Carrying out the Matsubara summations over Fermi frequencies and continuing an- 
alytically the result from the discrete points on the imaginary axis onto the complex 
plane, we find 



Appendix 




'k-kl +^k!-k 2 ~ U \ 



.0 , ,0 ,,oi 
'k-ki'^ki-ka'^kaJ 




(13) 
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The third-order self-energy Ti^\z) comprises two diagrams, a particle-hole one and 
a particle-particle one, 

4 3) W = 4 3 i(k,^) + 4 3 p(k,^) , (14) 
where 



ki,k 2 ,k 3 V k ~ kl ki-k 2 k 2 



Q(-<- k3 ,< 



Wki-ka ~ w k 3 J - ^k!-ka - w k 2 , 



(15) 



while the analogous expression for £p 3 p (k, z) is obtained from — £ P g(k, z) by chang- 
ing the sign of three a> k 's out of five, viz. w k _ kl , ^ kl _ k2 , and cf ki _ k3 , everything else 
being the same as in flT5]). The functions P and Q are defined as 

P(ei,e 2 ,e 3 ) = f( £l ) f(e 2 ) f(e 3 ) + f(-£i) f(-e 2 ) f(-e 3 ) , 
Q(e 1; e 2 ) = f( ei )f(e 2 )-f(- ei )f(-e 2 ) , 

f(e) = [e^ e + l)^ 1 is the Fermi function, with (3 = 1/ksT, and cj k = e k — //. 

The above equations are not suitable for numerical evaluation of £ k (a>±iO + ) 
and Ejf (&> ± i0 + ) because of two reasons. First, one would have to deal with four- 
dimensional and six-dimensional integrations in k space, respectively. Also, the 



denominators in expressions ( |T3|) and (IB) vanish along various closed contours in k 



space. The multiple integrals are therefore defined by their principal values and their 
evaluation requires a dense grid in k space, which makes the numerical procedures 
very time-consuming. 

The second- and third-order self-energies, given by Eqs. (|1^) - (0), are there- 
fore transformed here from the energy-momentum to the space-time representation, 
which allows an accurate and efficient evaluation by the fft algorithm. 

In order to establish transparent shorthand notation, we first define the Fourier 
transforms as the operators 



k 

duo 



N 

k 



oo 

itot 
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as well as the inverse transforms 



R 



dt e 



We also introduce shorthand for four- dimensional transforms 



T = -Fk->R Fw-*t an d T 1 = ^k^R Fw-^t 



— i 



The retarded/advanced self-energy in space-time representation is then defined as 



Sr/a(R,t) = .F [Ek^ ± z0" 



For the 2nd-order self-energy we start from the expression fll3|) for £ k (cj±«0 + ] 
and first evaluate the u) — *■ t Fourier transform. Using the relation 



uj ± i8 



which holds for any 5 > 0, we factorize the denominator in ([H|) into a product of 
three exponentials and find 



E^JR, t) = Ti U' 1 6(±t) S^\t) 



(16) 



(2) 

where (t) denotes the double convolution in momentum space, 
= E [A(k-k 1 ,t)A(k 1 -k 2 ,t)B*(k 2 ,t) 



ki,k 2 



B(k-k 1 ,t)B(k 1 -k 2 ,t) A*(k 2 ,t)] 



with functions A(k, t) and B(k, t) defined as 



A(k,f) 
B(k,t) 



e -K*f(_a;0) 



(17) 
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Recalling that the Fourier transform of a convolution can be expressed as a product 
of Fourier transforms of the integrands, we write A(k, t) and B(k, t) as 



A(M) 



[a(R,t)] 



(19) 



B(k,t) 



^ R [b(R,t)] 



(20) 



and disentangle the double convolution S£ 



(t) as 



S£> (t) = J^r [a 2 (R, t) b* (R, t) + b 2 (R, t) a* (R, t)] 



Inserted in expression (|TE|) this gives 



Sf/^R, t) = T« t^ 2 6(±t) [a 2 (R, i) b*(R, t) + b 2 (R, t) a*(R, t)] 



(21) 



In short, yields a convolution in k space, while .Fk-»R disentangles this convo- 
lution into a product. 

As regards the third-order self-energy, we proceed in the same way, starting 
from expression fll5|) for Spg(k, z). As before, the uo —>■ t Fourier transform fac- 
tories the ^-dependent term in the denominator, but the k summations still do 
not represent a convolution due to the presence of the unfactorized singular term 
l/(^ki-k 3 — ^3 — L<J ki-k2 + ^kJ- To factorize this term as well, we make use of the 



identity 




which holds for any e 7^ 0, and 




which holds for any 5 > 0, to represent 1/e as 




(22) 



s(f) 



sgn(t) e 
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It is clear that fl22"|) factorizes the above awkward term into a product of four expo- 
nentials. Altogether we obtain 



4 3 ^(R,t) = ±e(±t)t/ 3 / dt's(t')F^ n S^(k,t,t' 



A 3 ) , 



where the upper (lower) sign refers to the retarded (advanced) quantity, respectively, 
and 5pjj(k, t, t) denotes the triple k-space convolution 



AT3 



^ [A(k-k 1 ,t)B*(k 1 -k 2 ,t-t')A(k 2 ,t-t') 

ki,k 2 ,k 3 

+ B(k - k x , t) A*(ki - k 2 , t - t') B(k 2 , t - t')} 
x [A*(ki - k 3 , B(k 3 , - B*(k! - k 3 , A(k 3 , 0] • 



The particle-particle third-order term Ep P (k, z) is obtained from — Epg(k, z) by 
changing the sign of three cu^'s out of five, viz. a> k _ kl , w kl _ k2 , and ^; kl _ k3 . Since 
it is obvious from the above definitions (17 ) and (|18D of A(k, t) and B(k, t) that 



— c<J k implies A(k, t) «-> B*(k, i), one obtains Spp(k, £,£') from — Spj^k, t, t') 
by the appropriate replacements of those A's and B's that have k — ki, ki — k 2 , and 
k x — k 3 as their arguments: 



*Spp(k, £, t ) 



= "A73 E [B*(k-k 1 ,t)A(k 1 -k 2 ,t-t')A(k 2 ,t-t') 

ki,k 2 ,k 3 

+ A*(k - ki, t) B(ki - k 2 , t - t') B(k 2 , t - t')} 
x [B(kx - k 3 , B(k 3 , - A(kx - k 3 , A(k 3 , 0] • 

Expressing A(k, t) and B(k, t) in terms of their Fourier transforms, Eqs. ([19]) and 



(p20|) , as before, we decouple the above momentum-space convolutions and get 



E^R,*) = ±U 3 Q(±t) {[a(R,t)wi(R,t)-b(R,t)wJ(R,*)] 

- [a*(R,t)w 2 (R,t) +b*(R,t)w 3 (R,t)]} 



(23) 



where the functions Wj(R, t) are now given in the form of space-time convolutions 



w, 



/oo 
df^R-RV-OMR'O 
-°° 



(24) 
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with gj(R, t) and hj(R, t) denned as 



gi(R,t) = a(R,t)b*(R,t) , 

g 2 (R,t) = b 2 (R,t) , 

g 3 (R,t) = a 2 (R,t) , 

iu(R,t) = s(t)[gJ(R,t)-g 1 (R J t)] , 

h 2 (R,t) = s(t)[g 2 (R,t)-g 3 (R,t)] , 

h 3 (R,t) = h 2 (R,t) . 

For the 3rd order we thus need two more Fourier transforms than for the 2nd order, 
to decouple these additional space-time convolutions as well: 

w i (R,*)=^{^- 1 [g < (R,*)].^ 1 [h i (R,t)]} . (25) 
So, in order to obtain the nth-order self-energy E^(a;±iO + ) via its space-time 

n) 

1° 



representation E^(R, t) using the fft technique, we first calculate the quantities 



a(R,t) = ^[A^t^i^eM)^ , (26) 

k 

b(R,f) = ^ R [B(k,t)] = l^e'( k - R "^)f(-^) , (27) 

k 

which requires performing the Fourier transform jF k ^ R . Then we proceed to evaluate 
S^(R, t). For the 2nd order this amounts just to forming products of a(R, i)'s and 
b(R, i)'s, according to expression (plj). For the 3rd order we first form functions 
gi(R, t) and hj(R, t), again given by products of a(R, t)'s and b(R, £)'s, then perform 
a pair of Fourier transforms (JF _1 ,JF), to calculate functions w*(R, t) according to 
fl25|), and finally form £V(R, t) according to expression ([23|). Having thus found 
the self-energy in space-time representation, we apply once again to get 



The numerical problem involving space-time to momentum-frequency trans- 
formations (and vice versa) is solved by considering a finite lattice with periodic 
boundary conditions and discretizing the time axis. Here we take 256 x 256 lattice 
sites and 1024 time points. As regards the k «-> R transformation, the fft is an 
exact procedure, while the t <-> ui transformation involves approximating the con- 



13 



tinuous Laplace transform with the corresponding discrete fft. It would be much 
more time consuming to evaluate the space-time convolutions directly. 

The 4th-order self-energy (z) comprises twelve topologically nonequivalent 
diagrams, nine of which are numerically different. However, only four out of these 
nine can be treated entirely along the same lines as those of the 2nd and 3rd order, 
i.e. completely decoupled and evaluated by a series of Fourier transforms alone. The 
remaining five do not yield multiple space-time convolutions, but double integrals 
instead, which can be only partly decoupled by Fourier transforms. The 4th order 
thus introduces new numerical difficulties, but tractable ones. The fft can not do 
the whole job, but it does all steps but one. 

The functions a(R, t) and b(R, t), the building blocks of E^(R, t), have a 
clear physical interpretation as double-time correlation functions 

a(R,t) = <4(0)c R (t)> , 
b(R,t) = <c R (i)4(0)> , 

which can be shown by a straightforward evaluation of these thermal averages over 
the eigenstates of H , which gives and (|2"7|). This relates them to the unper- 
turbed Green's functions in space-time representation, viz. the retarded/advanced 
one, 

G° r/a (K,t) = T iG(±t) <[4(0),c R (t)] + > 
= =FiQ(±*)[b(R,t) + a(R,t)] , 

and the causal one, 

G°(R,t) = -i <T{c R (t) 4(0)} > 

= -i [6(t) b(R, t) - B(-t) a(R, *)] . 

It is clear that b(R, t) is the probability amplitude for an electron to be created at 
R = and t = 0, to propagate to the site R which is reached at a later time t, and to 
be destructed at (R, t). This newly created electron is "composed of" the unoccupied 
(i.e., still available) k states, as indicated by the amplitudes f(— u;£) = 1— f(u;£) in the 
defining relation (57]). In the same way a(R, t) effectively describes the propagation 
of a hole created at (R, t < 0) and destructed at (R = 0, t — 0). 
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Figure captions 



Fig. 1. Diagrammatic representation of second- and third-order self-energy contri- 
butions. 

Fig. 2. Real part of E^(u;) along the r — > X, X — > Z, and Z — > T cuts through the 
Brillouin zone. ReEp 3 ^, — ReEp P , and ReE® = ReEp 3 ^ + ReE p 3 p are represented 
by dashed, dashed-dotted, and full lines, respectively. 

Fig. 3. Imaginary part of E k (cu) along the r — > X, X — ■> Z, and Z — > T cuts through 
the Brillouin zone. ImE p 3 ^, — ImE p 3 p , and Im = Im Ep 3 -^ + Im Ep 3 p are repre- 
sented by dashed, dashed-dotted, and full lines, respectively. 

Fig. 4. Imaginary parts of S ( JV) an d £ ( k M alon S the r X, X Z, and Z -> T 
cuts through the Brillouin zone. ImS( 2 ', ImS®, and the sum ImS^ + ImE® are 
represented by dashed, dashed-dotted, and full lines, respectively. 

Fig. 5. Imaginary parts of E^(cj) and E^(a>) at T, X, Z, and M points close to 
uo = E F . ImS' 2 ', ImS®, and the sum Im S* 2 ' + Im S* 3 ' are represented by dashed, 
dashed-dotted, and full lines, respectively. 

Fig. 6. Single-particle spectral function Ak(co>) along the r — > X, X — > Z, and Z — > T 
cuts through the Brillouin zone. 

Fig. 7. Quasiparticle dispersion (circles) derived from the maxima of the spectral 
function Ak(o>) and the unrenormalized dispersion (full line). 

Fig. 8. Renormalized (full line) and unrenormalized (dashed line) density of states. 
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Figure V. Diagrammatic representation of second- and third-order self-energy con- 
tributions. 
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Figure 2: Real part of S k (a;) along the r — > X, X — > Z, and Z — > T cuts through the 
Brillouin zone. Re£ pp: , — Re£ pp , and ReS^ 3 ) = Re£ pp: + Re£ pp are represented 
by dashed, dashed-dotted, and full lines, respectively. 
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Figure 3: Imaginary part of E k (a;) along the r — > X, X — > Z, and Z — > T cuts 
through the Brillouin zone. ImSp 3 ^, — ImEp 3 p , and ImS® = ImEp 3 ^ + ImSpp are 
represented by dashed, dashed-dotted, and full lines, respectively. 
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Figure 4: Imaginary parts of X k (u) and E k \oj) along the r — > X, X — > Z, 
and Z — > T cuts through the Brillouin zone. ImS( 2 \ ImE®, and the sum 
ImS' 2 ' +ImS^ are represented by dashed, dashed-dotted, and full lines, respec- 
tively. 
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Figure 5: Imaginary parts of £ k (uj) and £ k (uj) at T, X, Z, and M points close to 
uj = Ep. ImE' 2 \ ImE®, and the sum Im E' 2 ' + Im E' 3 ' are represented by dashed, 
dashed-dotted, and full lines, respectively. 
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Figure 6: Single-particle spectral function Ak(cj) along the r — > X, X — > Z, and 
Z — > T cuts through the Brillouin zone. 
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Figure 7: Quasiparticle dispersion (circles) derived from the maxima of the spectral 
function Aj c (o;) and the unrenormalized dispersion (full line). 




-8 -4 4 8 

UJ I t 



Figure 8: Renormalized (full line) and unrenormalized (dashed line) density of states. 
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